///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/wfn/configuration.cpp
///
/// \brief   implementation of members for the Configuration and ConfigDiff classes
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include<algorithm>
#include<numeric>

#include<formic/wfn/configuration.h>
#include<formic/random/random.h>
#include<formic/input/base.h>

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief default constructor
///
///////////////////////////////////////////////////////////////////////////////////////////////////
formic::ConfigDiff::ConfigDiff()
  : _newly_empty(),
    _newly_full()
{}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Constructs the configuration difference object
///
/// \param[in]      userinp  the user's input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////
formic::ConfigDiff::ConfigDiff(const formic::InputBase & userinp)
  : _newly_empty(),
    _newly_full()
{}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  constructor from two configurations
///
/// \param[in]      new_config     the new configuration
/// \param[in]      old_config     the old configuration
///
///////////////////////////////////////////////////////////////////////////////////////////////////
formic::ConfigDiff::ConfigDiff(const formic::Configuration & new_config, const formic::Configuration & old_config)
  : _newly_empty(),
    _newly_full()
{

  // ensure the fock vectors are the same lengths
  if ( new_config.fock_config().size() != old_config.fock_config().size() )
    throw formic::Exception("fock vector size mismatch in ConfgDiff constructor");

  // process alpha and beta sites separately
  for (int b = 0; b < 2; b++) {

    // check for sites for which the occupation changed
    for (int i = b; i < new_config.fock_config().size(); i+=2)
      if (new_config[i] == 0 && old_config[i] == 1)
        _newly_empty.push_back(i);
      else if (new_config[i] == 1 && old_config[i] == 0)
        _newly_full.push_back(i);

    // ensure there were the same number of particles created as destroyed
    if ( _newly_empty.size() != _newly_full.size() )
      throw formic::Exception("uneven number of newly empty/full sites for %s particles in ConfgDiff constructor")
            % ( b == 0 ? "alpha" : "beta" );

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  return a compound index for a pair of sites, one occupied and the other empty
///
/// \param[in]     occ_site     the index of the occupied site (in alpha and beta indexing)
/// \param[in]     unocc_site   the index of the unoccupied site (in alpha and beta indexing)
///
///////////////////////////////////////////////////////////////////////////////////////////////////
int formic::Configuration::occ_unocc_cmpd(int occ_site, int unocc_site) const {

  assert( _cps == 2 );

  const int b = occ_site % 2;

  assert( occ_site % 2 == unocc_site % 2 );

  const std::vector<int> & site_to_occ   = ( b == 0 ? _a_site_to_occ   : _b_site_to_occ   );
  const std::vector<int> & site_to_unocc = ( b == 0 ? _a_site_to_unocc : _b_site_to_unocc );
  const int np = ( b == 0 ? _nap : _nbp );
  const int nu = ( b == 0 ? _nau : _nbu );

  occ_site /= 2;
  unocc_site /= 2;

  assert( site_to_occ[occ_site] >= 0 );
  assert( site_to_occ[occ_site] < np );

  int retval;

  if ( occ_site == unocc_site ) {

    retval = site_to_occ[occ_site];

  } else {

    assert( site_to_unocc[unocc_site] >= 0 );
    assert( site_to_unocc[unocc_site] < nu );

    retval = np + nu * site_to_occ[occ_site] + site_to_unocc[unocc_site];

  }

  assert( retval >= 0 );
  assert( retval < np*(nu+1) );

  return retval;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief proposes a move from one configuration to another
///
/// \param[in]     userinp      the user's input options
/// \param[in]     gen          a random number generator
/// \param[in]     config       the configuration to move from
/// \param[out]    config_diff  on exit, specifies the difference between the configurations
///
/// \return the proposal density ratio, which is the ratio of the probability of proposing the
///         current configuration if we are at the new configuration to the probability
///         of proposing the new configuration if we are at the current configuration.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class GEN> double formic::Configuration::propose_move(const formic::InputBase & userinp,
                                                                GEN & gen,
                                                                const formic::Configuration & config,
                                                                formic::ConfigDiff & config_diff) {

  // get number of atomic orbitals
  const int nao = userinp.get<int>("nao");

  // check that this is a 2-configuration-per-site configuration with an even number of sites
  if ( config.cps() != 2 || config.fock_config().size() != 2*nao )
    throw formic::Exception("formic::Configuration::propose_move requires a fermionic configuration object");

  // propose a new fock space configuration by moving one alpha or one beta electron
  const int pos1 = gen() % userinp.get<int>("nsites");
  const int pos1_occ = config.fock_config().at(pos1);
  assert(pos1_occ == 0 || pos1_occ == 1);
  const int b = pos1 % 2;
  const int n_a_part = ( b == 0 ? userinp.get<int>("nalpha") : userinp.get<int>("nbeta") );
  const int x = gen() % ( pos1_occ == 1 ? nao - n_a_part : n_a_part );
  int c = 0;
  int j = 0;
  for ( ; j < nao; j++)
    if (config.fock_config().at(2*j+b) != pos1_occ)
      if (c++ == x)
        break;
  const int pos2 = 2*j+b;
  const int pos2_occ = config.fock_config().at(pos2);
  assert(c == x+1);
  assert(pos2_occ == 0 || pos2_occ == 1);
  assert(pos1 != pos2);
  assert(pos1_occ != pos2_occ);
  assert(pos1 % 2 == pos2 % 2);

  // define the difference
  config_diff.newly_empty().assign(1, (pos1_occ == 1 ? pos1 : pos2));
  config_diff.newly_full().assign(1, (pos1_occ == 0 ? pos1 : pos2));

  // for now, the proposal density is symmetric, so return one
  return 1.0;

}

template double formic::Configuration::propose_move(const formic::InputBase & userinp,
                                                    formic::LaggedFibonacci<> & gen,
                                                    const formic::Configuration & config,
                                                    formic::ConfigDiff & config_diff);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Constructs a fermionic configuration object from the user's input options
///
/// \param[in]      userinp  the user's input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////
formic::Configuration::Configuration(const formic::InputBase & userinp)
  : _ns(2*userinp.get<int>("nao")),
    _to( userinp.get<int>("nalpha") + userinp.get<int>("nbeta") ),
    _cps(2),
    _finished(false),
    _occ_vec(userinp.get<int>("nsites")),
    _occurances(2, 0),
    _nas(userinp.get<int>("nao")),
    _nbs(userinp.get<int>("nao")),
    _nap(userinp.get<int>("nalpha")),
    _nbp(userinp.get<int>("nbeta")),
    _nau( _nas - _nap ),
    _nbu( _nbs - _nbp ),
    _a_site_to_occ(_nas, -1),
    _b_site_to_occ(_nbs, -1),
    _a_site_to_unocc(_nas, -1),
    _b_site_to_unocc(_nbs, -1),
    _a_occ_to_site(_nap, -1),
    _b_occ_to_site(_nbp, -1),
    _a_unocc_to_site(_nau, -1),
    _b_unocc_to_site(_nbu, -1)
{

  if ( _nap != _nbp )
    throw formic::Exception("configuration constructor currently assumes equal number of alpha and beta");

  // set some constants
  const int nsites = userinp.get<int>("nsites");
  std::vector<int> ab_n_particles(2, 0);
  ab_n_particles.at(0) = userinp.get<int>("nalpha");
  ab_n_particles.at(1) = userinp.get<int>("nbeta");

  // if supplied, start in the user's chosen fock space configuration
  if (userinp.get<std::vector<int> >("init_config").size() == nsites) {

    _occ_vec = userinp.get<std::vector<int> >("init_config");

  // otherwise, generate a random fock space configuration
  } else {

    // remove all particles from the configuration
    _occ_vec.assign(nsites, 0);

    // place alpha (b == 0) and beta (b == 1) particles
    for (int b = 0; b <= 1; b++) {
      for (int i = 0; i < ab_n_particles.at(b); i++) {
        int x = formic::global_lcg() % (nsites/2 - i);
        int c = 0;
        int j = 0;
        for ( ; j < nsites/2; j++)
          if (_occ_vec.at(2*j+b) == 0)
            if (c++ == x)
              break;
        assert(c == x+1);
        assert(_occ_vec.at(2*j+b) == 0);
        _occ_vec.at(2*j+b) = 1;
      }
    }

  }

  // initialize the occurances vector
  for (std::vector<int>::const_iterator it = _occ_vec.begin(); it != _occ_vec.end(); it++)
    _occurances.at(*it) += 1;

  // initialize occupation/unoccupation/site mappings
  this->init_occ_unocc_site();

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  assignment operator for the Configuration class
///
/// \param[in]      other    the configuration to copy from
///
/// \return  a reference to the object
///
///////////////////////////////////////////////////////////////////////////////////////////////////
formic::Configuration & formic::Configuration::operator=(const formic::Configuration & other) {

  if (&other == this)
    return *this;

  _ns = other._ns;
  _to = other._to;
  _cps = other._cps;
  _finished = other._finished;

  _nas = other._nas;
  _nbs = other._nbs;
  _nap = other._nap;
  _nbp = other._nbp;
  _nau = other._nau;
  _nbu = other._nbu;

  if ( _occ_vec.size() != other._occ_vec.size() ) _occ_vec.resize(other._occ_vec.size());
  if ( _occurances.size() != other._occurances.size() ) _occurances.resize(other._occurances.size());
  if ( _a_site_to_occ.size() != other._a_site_to_occ.size() ) _a_site_to_occ.resize(other._a_site_to_occ.size());
  if ( _b_site_to_occ.size() != other._b_site_to_occ.size() ) _b_site_to_occ.resize(other._b_site_to_occ.size());
  if ( _a_site_to_unocc.size() != other._a_site_to_unocc.size() ) _a_site_to_unocc.resize(other._a_site_to_unocc.size());
  if ( _b_site_to_unocc.size() != other._b_site_to_unocc.size() ) _b_site_to_unocc.resize(other._b_site_to_unocc.size());
  if ( _a_occ_to_site.size() != other._a_occ_to_site.size() ) _a_occ_to_site.resize(other._a_occ_to_site.size());
  if ( _b_occ_to_site.size() != other._b_occ_to_site.size() ) _b_occ_to_site.resize(other._b_occ_to_site.size());
  if ( _a_unocc_to_site.size() != other._a_unocc_to_site.size() ) _a_unocc_to_site.resize(other._a_unocc_to_site.size());
  if ( _b_unocc_to_site.size() != other._b_unocc_to_site.size() ) _b_unocc_to_site.resize(other._b_unocc_to_site.size());

  std::copy(other._occ_vec.begin(), other._occ_vec.end(), _occ_vec.begin());
  std::copy(other._occurances.begin(), other._occurances.end(), _occurances.begin());
  std::copy(other._a_site_to_occ.begin(), other._a_site_to_occ.end(), _a_site_to_occ.begin());
  std::copy(other._b_site_to_occ.begin(), other._b_site_to_occ.end(), _b_site_to_occ.begin());
  std::copy(other._a_site_to_unocc.begin(), other._a_site_to_unocc.end(), _a_site_to_unocc.begin());
  std::copy(other._b_site_to_unocc.begin(), other._b_site_to_unocc.end(), _b_site_to_unocc.begin());
  std::copy(other._a_occ_to_site.begin(), other._a_occ_to_site.end(), _a_occ_to_site.begin());
  std::copy(other._b_occ_to_site.begin(), other._b_occ_to_site.end(), _b_occ_to_site.begin());
  std::copy(other._a_unocc_to_site.begin(), other._a_unocc_to_site.end(), _a_unocc_to_site.begin());
  std::copy(other._b_unocc_to_site.begin(), other._b_unocc_to_site.end(), _b_unocc_to_site.begin());

  return *this;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Initializes the alpha (b = 0) or beta (b = 1) occ <--> site and unocc <--> site
///         mappings for the stored configuration
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::Configuration::init_occ_unocc_site(const int b) {

  std::vector<int> & site_to_occ   = ( b == 0 ? _a_site_to_occ   : _b_site_to_occ   );
  std::vector<int> & site_to_unocc = ( b == 0 ? _a_site_to_unocc : _b_site_to_unocc );
  std::vector<int> & occ_to_site   = ( b == 0 ? _a_occ_to_site   : _b_occ_to_site   );
  std::vector<int> & unocc_to_site = ( b == 0 ? _a_unocc_to_site : _b_unocc_to_site );
  const int ns = ( b == 0 ? _nas : _nbs );
  const int np = ( b == 0 ? _nap : _nbp );
  const int nu = ( b == 0 ? _nau : _nbu );

  site_to_occ.resize(ns);
  site_to_unocc.resize(ns);
  occ_to_site.resize(np);
  unocc_to_site.resize(nu);

  assert( _occ_vec.size() == 2*ns );

  int pc = 0; // particle counter
  int hc = 0; // hole counter

  int i = b;
  for ( ; i < _occ_vec.size(); i+=2) {

    if ( _occ_vec[i] == 1 ) {

      site_to_unocc[i/2] = -1;
      site_to_occ[i/2] = pc;
      occ_to_site[pc] = i/2;
      pc++;

    } else {

      site_to_occ[i/2] = -1;
      site_to_unocc[i/2] = hc;
      unocc_to_site[hc] = i/2;
      hc++;

    }

  }

  assert( i/2 == ns );
  assert(  pc == np );
  assert(  hc == nu );

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Initializes the occ <--> site and unocc <--> site mappings for the stored configuration
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::Configuration::init_occ_unocc_site() {

  this->init_occ_unocc_site(0);
  this->init_occ_unocc_site(1);

//  for (int b = 0; b < 2; b++) {
//
//    std::vector<int> & site_to_occ   = ( b == 0 ? _a_site_to_occ   : _b_site_to_occ   );
//    std::vector<int> & site_to_unocc = ( b == 0 ? _a_site_to_unocc : _b_site_to_unocc );
//    std::vector<int> & occ_to_site   = ( b == 0 ? _a_occ_to_site   : _b_occ_to_site   );
//    std::vector<int> & unocc_to_site = ( b == 0 ? _a_unocc_to_site : _b_unocc_to_site );
//    const int ns = ( b == 0 ? _nas : _nbs );
//    const int np = ( b == 0 ? _nap : _nbp );
//    const int nu = ( b == 0 ? _nau : _nbu );
//
//    site_to_occ.resize(ns);
//    site_to_unocc.resize(ns);
//    occ_to_site.resize(np);
//    unocc_to_site.resize(nu);
//
//    assert( _occ_vec.size() == 2*ns );
//
//    int pc = 0; // particle counter
//    int hc = 0; // hole counter
//
//    int i = b;
//    for ( ; i < _occ_vec.size(); i+=2) {
//
//      if ( _occ_vec[i] == 1 ) {
//
//        site_to_unocc[i/2] = -1;
//        site_to_occ[i/2] = pc;
//        occ_to_site[pc] = i/2;
//        pc++;
//
//      } else {
//
//        site_to_occ[i/2] = -1;
//        site_to_unocc[i/2] = hc;
//        unocc_to_site[hc] = i/2;
//        hc++;
//
//      }
//
//    }
//
//    assert( i/2 == ns );
//    assert(  pc == np );
//    assert(  hc == nu );
//
//  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Updates the object to incorperate the supplied difference.
///
/// \param[in]     diff  specifies the difference between the configuration to update to
///                      and the current configuration
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::Configuration::update(const formic::ConfigDiff & diff) {

  assert( diff.newly_full().size() == diff.newly_empty().size() );

  for (int i = 0; i < diff.newly_full().size(); i++) {

    int new_occ = diff.newly_full().at(i);
    int old_occ = diff.newly_empty().at(i);

    assert( new_occ % 2 == old_occ % 2 );

    assert( _occ_vec.at(new_occ) == 0 );
    assert( _occ_vec.at(old_occ) == 1 );

    _occ_vec.at(new_occ) = 1;
    _occ_vec.at(old_occ) = 0;

    const int b = new_occ % 2;

    std::vector<int> & site_to_occ   = ( b == 0 ? _a_site_to_occ   : _b_site_to_occ   );
    std::vector<int> & site_to_unocc = ( b == 0 ? _a_site_to_unocc : _b_site_to_unocc );
    std::vector<int> & occ_to_site   = ( b == 0 ? _a_occ_to_site   : _b_occ_to_site   );
    std::vector<int> & unocc_to_site = ( b == 0 ? _a_unocc_to_site : _b_unocc_to_site );
    const int ns = ( b == 0 ? _nas : _nbs );
    const int np = ( b == 0 ? _nap : _nbp );
    const int nu = ( b == 0 ? _nau : _nbu );

    new_occ /= 2;
    old_occ /= 2;

    assert( site_to_occ[old_occ] >= 0 );
    assert( site_to_occ[old_occ] < np );
    assert( occ_to_site[ site_to_occ[old_occ] ] == old_occ );
    occ_to_site[ site_to_occ[old_occ] ] = new_occ;
    site_to_occ[new_occ] = site_to_occ[old_occ];
    site_to_occ[old_occ] = -1;

    assert( site_to_unocc[new_occ] >= 0 );
    assert( site_to_unocc[new_occ] < nu );
    assert( unocc_to_site[ site_to_unocc[new_occ] ] == new_occ );
    unocc_to_site[ site_to_unocc[new_occ] ] = old_occ;
    site_to_unocc[old_occ] = site_to_unocc[new_occ];
    site_to_unocc[new_occ] = -1;

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief checks that using the supplied difference object to update the supplied configuration
///        is reasonable
///
/// \param[in]      diff          object describing a change to the configuration
/// \param[in]      config        the configuration
/// \param[in]      calling_func  the name of the function in which this check is performed
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::check_config_update(const formic::ConfigDiff & diff,
                                 const formic::Configuration & config,
                                 const std::string & calling_func) {

  // check newly full sites
  for (int i = 0; i < diff.newly_full().size(); i++) {
    if ( config[diff.newly_full().at(i)] != 0 )
      throw formic::Exception("check_config_update failed: cannot move an electron into a site that is already occupied in %s") % calling_func;
    for (int j = 0; j < i; j++)
      if ( diff.newly_full().at(i) == diff.newly_full().at(j) )
        throw formic::Exception("check_config_update failed: cannot move two electrons into the same site in %s") % calling_func;
  }

  // check newly empty sites
  for (int i = 0; i < diff.newly_empty().size(); i++) {
    if ( config[diff.newly_empty().at(i)] == 0 )
      throw formic::Exception("check_config_update failed: cannot move an electron out of a site that is empty in %s") % calling_func;
    for (int j = 0; j < i; j++)
      if ( diff.newly_empty().at(i) == diff.newly_empty().at(j) )
        throw formic::Exception("check_config_update failed: cannot move two electrons out of the same site in %s") % calling_func;
  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Constructs the a configuration object with arbitrary configurations per site
///
/// \param[in]      ns       the number of sites
/// \param[in]      to       the total lattice occupation
/// \param[in]      cps      the number of configurations on a site
///
///////////////////////////////////////////////////////////////////////////////////////////////////
formic::Configuration::Configuration(const int ns, const int to, const int cps)
  : _ns(ns),
    _to(to),
    _cps(cps),
    _occ_vec(ns, 0),
    _occurances(cps, 0)
{

//  // get the number of configurations, using the occurances vector as workspace
//  _nc = formic::n_integer_solutions(ns, to, cps-1, &_occurances.at(0));
//
//  // reset the occurances vector to zero
//  std::fill(_occurances.begin(), _occurances.end(), 0);

  // reset to the first valid configuration
  this->reset();

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Resets the configuration
///        (currently only set up for the arbitrary configs per site portion of the object)
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::Configuration::reset() {

  // reset the finished flag
  _finished = false;

  // starting with the largest occupation, set each occupation's occurance to its smallest
  // possible value
  {
    int rem_occ = _to;
    int rem_sites = _ns;
    for (int i = _cps - 1; i > 0; i--) {
      _occurances[i] = 0;
      while( rem_occ - i * _occurances[i] - (i-1) * (rem_sites - _occurances[i]) > 0 )
        _occurances[i]++;
      rem_sites -= _occurances[i];
      rem_occ -= i * _occurances[i];
    }
    assert( rem_occ == 0 );
    assert( rem_sites >= 0 );
  }

  // sort the occupations into the occupation number vector in descending order
  {
    int pos = 0;
    for (int i = _cps - 1; i > 0; i--)
      for (int j = 0; j < _occurances[i]; j++)
        _occ_vec[pos++] = i;
    for ( ; pos < _occ_vec.size(); pos++)
      _occ_vec[pos] = 0;
  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Increments to the next valid configuration
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::Configuration::operator++() {

  // loop over the different particle ranks
  bool finished_shifting = true;
  int r = 1;
  while ( r < _cps ) {

    // shift right
    finished_shifting = false;
    bool flag = false;
    int accumulation = 0;
    for (int i = _occ_vec.size() - 1; i >= 0; i--) {
      
      if (_occ_vec[i] == r)
        accumulation++;

      if ( accumulation == _occurances[r] && !flag ) {
        finished_shifting = true;
        break;
      }

      if (_occ_vec[i] < r)
        flag = true;

      if (_occ_vec[i] == r && flag) {
        for (int k = 0; k < i; k++)
          if (_occ_vec[k] < r)
            _occ_vec[k] = -1;
        _occ_vec[i] = -1;
        //accumulation--;
        int pos = i+1;
        for ( ; accumulation > 0; pos++) {
          assert( pos < _occ_vec.size() );
          if (_occ_vec[pos] <= r) {
            _occ_vec[pos] = r;
            accumulation--;
          }
        }
        assert( pos <= _occ_vec.size() );
        for ( ; pos < _occ_vec.size(); pos++)
          if (_occ_vec[pos] <= r)
            _occ_vec[pos] = -1;
        break;
      }

    }

    // if we are not finished shifting at this particle rank, then the increment is complete
    if (!finished_shifting) break;

    // otherwise, we re-distribute the next particle rank in the next iteration of this loop
    r++;

  }

  // if we were are finished shifting, move to the next distribution of sites among particle ranks
  if (finished_shifting) {

    // set the occupation vector to -1 flags
    for (int i = 0; i < _occ_vec.size(); i++)
      _occ_vec[i] = -1;

    // get the remaining occupation
    int rem_occ = _to;
    int rem_sites = _ns;
    for (int s = _cps-1; s > 1; s--) {
      rem_sites -= _occurances[s];
      rem_occ -= s * _occurances[s];
    }

    // redistribute
    int s = 2;
    for ( ; s < _cps; s++) {
      // try adding another occupancy of this rank
      if ( rem_occ - s >= 0 && rem_sites - 1 >= 0 ) {
        _occurances[s]++;
        rem_sites--;
        rem_occ -= s;
        break;
      }
      // if that doesn't work, go up to the next particle rank
      rem_sites += _occurances[s];
      rem_occ += s * _occurances[s];
    }
    //std::cout << s << std::endl;
    if (s == _cps) {
      _finished = true;
      return;
    }
    for (s = s - 1; s > 0; s--) {
      _occurances[s] = 0;
      while( rem_occ - s * _occurances[s] - (s-1) * (rem_sites - _occurances[s]) > 0 )
        _occurances[s]++;
      rem_sites -= _occurances[s];
      rem_occ -= s * _occurances[s];
    }
    assert( rem_occ == 0 );
    assert( rem_sites >= 0 );

    // make sure that r is maximum so that all particle rank's arrangements will be reset
    assert( r == _cps );

  }

  //printf("  %10s    ", " ");
  //for (int i = 0; i < _occ_vec.size(); i++)
  //  printf(" %2i", _occ_vec.at(i));
  //printf("    %2i    %2i", r, _occurances[1]);
  //std::cout << std::endl;

  // reset the lower particle ranks now that we have shifted or moved to a new distribution
  {
    int accumulation = 0;
    for (int pos = 0; pos < _occ_vec.size(); pos++) {
      if (_occ_vec[pos] == -1) {
        while ( accumulation == 0 && r > 0 ) {
          r--;
          accumulation = _occurances[r];
        }
        _occ_vec[pos] = r;
        accumulation--;
      }
    }
    //assert( r == 0 );
  }

  //printf("  %10s    ", " ");
  //for (int i = 0; i < _occ_vec.size(); i++)
  //  printf(" %2i", _occ_vec.at(i));
  //std::cout << std::endl;

  // check that the total occupation remains correct
  assert( std::accumulate(_occ_vec.begin(), _occ_vec.end(), 0) == _to );

  //assert( 1 == 0 );

}

/// \brief  function to populate the alpha or beta configuration based on the compound index x
void formic::Configuration::populate_alpha_or_beta(const int b, int x) {

  assert( b == 0 || b == 1 );
  assert( _cps == 2 );
  assert( _occ_vec.size() % 2 == 0 );
  const int no = ( b == 0 ? _nas : _nbs );
  assert( _occ_vec.size() == 2*no );
  int & np = ( b == 0 ? _nap : _nbp );
  np = 0;
  int * fin = &_occ_vec[b] + 2*no;
  for (int * occ = &_occ_vec[b]; occ != fin; occ += 2) {
    _occurances[ *occ ] -= 1;
    const int new_occ = x - 2 * ( x / 2 ); //x % 2;
    np += new_occ;
    *occ = new_occ;
    _occurances[ *occ ] += 1;
    x /= 2;
  }
  int & nu = ( b == 0 ? _nau : _nbu );
  nu = no - np;
  assert( nu >=  0 );
  assert( nu <= no );
  assert( np >=  0 );
  assert( np <= no );
  this->init_occ_unocc_site(b);

}

/// \brief  returns the irreducible representation of the configuration
int formic::Configuration::get_irrep(const int nirr, const std::vector<int> & orb_irrep, const std::vector<int> & dir_prod_table) const {
  int retval = 0;
  for (int i = 0; i < _occ_vec.size(); i++)
    retval = ( _occ_vec[i] == 0 ? retval : dir_prod_table.at( orb_irrep.at(i/2) + nirr * retval ) );
  return retval;
}
